Relatório Parcial

Alinhamento de Sequências de DNA

Alinhamento de sequências de DNA consiste na comparação entre duas sequências a fim de, por exemplo, tentar relacionar estes dados a funções biológicas ou estruturas presentes nos organismos. Para este estudo, foram utilizados 3 algoritmos: Smith-Waterman, Busca Local e Busca Exaustiva. Para a busca local foi feita com base no conceito de aleatoriedade, artifício muito válido quando utilizado em conjunto com outros. Com relação a busca exaustiva, foram feitas também algumas otimizações.

Este relatório parcial tem como objetivo analisar os tempos de execução destes algoritmos para diferentes tamanhos de entrada, além dos possíveis pontos de lentidão. Os códigos e arquivos para execução estão presentes no repositório. Foram utilizados 2 testes para cada tamanho, totalizando em 200 testes.

Uma breve explicação sobre cada algoritmo:

Algoritmo para alinhamento de sequências de DNA exaustiva

Sendo i o índice do início da substring A, j o índice do início da substring B, k o tamanho das substrings, n o tamanho da sequência A, m o tamanho da sequência B, subA a substring A e subB a substring B, utiliza-se o pseudo-código abaixo para fazer os alinhamentos:

  1. k = tamanho da menor sequência

  2. Enquanto k >= 2 e k > tamanho mínimo de substring que pode atingir o maior score:

    1. Enquanto i < n-k:

      1. subA = substring de tamanho k a partir do indice i

      2. Enquanto j < m-k:

        1. subB = substring de tamanho k a partir do indice i

        2. Calcular score entre subA e subB

O código no estado atual não percorre todas as possibilidades como um algoritmo exaustivo. Isso ocorre, pois foram efetuadas algumas modificações. A primeira verifica substrings repetidas, ou seja, para cada tamanho k, as substrings verificadas para cada sequência são armazenadas para comparação com as novas substrings geradas. Caso uma esteja previamente armazenada, os loops subsequentes não serão efetuados. Além disso, a partir do score máximo, é calculado o tamanho de string mínimo para que seja possível atingir um score igual ou maior, ou seja, como o valor para um match é de 2 pontos, caso o score máximo em certo ponto do algoritmo seja 16, então o tamanho mínimo será 8.

Algoritmo de Smith Waterman

O algoritmo tem a finalidade de fazer alinhamento local de duas sequências de DNA. Ele utiliza como entrada as duas sequências e o tamanho delas. Sua saída é um arquivo .txt contendo a matriz com os valores, o maior valor encontrado e ambas as sequências alinhadas. A abordagem procura encontrar subsequências que maximizam os valores de pareamento. Para calcular a matriz de valores e obter o valor máximo, segue-se o seguinte algoritimo:

  1. Inicializar H[i,0]=0, 0≤i≤n

  2. Inicializar H[0,j]=0, 0≤j≤m

  3. Para cada 1≤i≤n e 1≤j≤m:

    1. Calcular diagonal = H[i-1,j-1] + w(a[i],b[j]),onde w(a[i],b[j])=2 se houve match, w(a[i],b[j])= -1 se houve mismatch ou gap

    2. Calcular deleção = H[i-1,j] - 1

    3. Calcular inserção = H[i,j-1] - 1

    4. Calcular H[i,j]=máximo (0, diagonal, deleção, inserção)

9 Retornar o máximo de H

Durante o calculo do máximo entre diagonal, deleção e inserção no passo 7, a origem do valor é guardada no H[i,j], ou seja, caso o máximo seja o da diagonal, então a origem guardada é um código referente a diagonal, caso seja de deleção, um código referente a deleção e caso seja de inserção, um código referente a inserção. Para obter as sequências alinhadas pode-se utilizar essas origens dos valores para gerar a sequência da frente para trás e invertê-la no final.

Algoritmo para alinhamento de sequências de DNA com aleatoriedade

Sendo i o índice do início da substring A, j o índice do início da substring B, n o tamanho da substring A, m o tamanho da substring B, k o tamanho das substrings e p a quantidade de substrings geradas em B, utiliza-se o pseudo-código abaixo para fazer os alinhamentos:

  1. Por 10000 vezes: (valor arbitrário para agilizar processo de análise do algoritmo)

    1. Se n>m, gerar k de 0 até m

    2. Senão, gerar k de 0 até n

    3. Gerar p subsequencias sa=a[i,i+1,...,i+k] de a, com tamanho k, 0<=i<=n

    4. Gerar um número aleatório inteiro positivo p entre 0 e m

    5. Por p vezes:

      1. Gerar uma subsequencia sb=b[j,j+1,...,j+k] de b, de tamanho aleatório k, 1<=k<=m, e 0<=j<=m

      2. Calcular os scores de cada par (sa,sb) com os pesos wmat, wmis e wgap

9 Devolver o score máximo m entre os scores do passo (4) e as subsequencias associadas a ele

Alguns tratamentos foram efetuados para que o programa funcione corretamente. Primeiramente verifica-se qual a menor sequência entre as duas que serviram de entrada para o cálculo do valor k, evitando que este seja maior que o tamanho de uma das sequências originais. Outro tratamento é que os valores de j devem ser entre 0 e m-k. Os valores de j, k e p são gerados aleatóriamente. Além disso, os valores de wmat, wmis, wgap são os mesmos utilizados no outro projeto (https://github.com/lucaskf1996/SmithWaterman) para que uma comparação possa ser feita eventualmente. Outras obsevações a serem feitas é que o cálculo do score é feito de caractere em caractere, podendo ter um resultado negativo e que i foi mantido como 0, mas que pode ser facilmente modificado para adicionar mais aleatoriedade à solução.

Análise de tempo de execução

Os tempos de execução foram extraídos utilizando o código python (teste.py) presente no repositório. Estes foram escritos em arquivos .txt e tratados no código abaixo.

Preparando as variáveis para os gráficos:

Com os dados tratados, foram feitos gráficos 3D através da biblioteca plotly:

Smith-Waterman

No gráfico é possível observar certos picos no tempo de execução. Para tentar esclarecer o motivo para tal efeito, foram feitas novas baterias de testes utilizando o mesmo código com os mesmos testes.

Como é possível observar, os picos de tempo de execução não estão presentes nos 3 gráficos ao mesmo tempo. Isto é um indicativo de que possívelmente não é um problema de código em si. Uma possibilidade é que isso seja um problema com o tempo limite do uso da thread alocada para o processo, ou seja, o tempo de uso da thread acaba durante a execução do programa e este é então alocado a outra thread, consequentemente, aumentando o tempo de execução.

Algoritmo Local

Algoritmo Exaustivo

Observando os gráficos fica claro que todos estão dependentes dos tamanhos das entradas. Entre os 3 algoritmos, o Smith-Waterman foi o com o menor tempo de execução, além dos melhores scores, algo que teóricamente seria impossível, considerando que um dos outros algoritmos é um exaustivo, ou seja, o score deveria ser no máximo igual. O fato dos scores serem melhores pode ser explicado por pelo fato de que o algoritmo ajusta as strings resultantes caso haja uma inserção ou deleção, efeito não reproduzido nos outros dois. Quanto ao tempo de execução, todos os algoritmos dependem dos tamanhos das entradas, sendo a exaustiva a que mais escala, porém esta, quando comparada ao algoritmo local, possui um score muito melhor.

Abaixo estão representadas os tempos de execução para entradas de de mesmo tamanho n e m:

Analisando os códigos é possível identificar que tanto os algoritmos local e Smith-Waterman são ambos O(n²) por conta dos dois loops aninhados. Para o Smith-Waterman, há dois loops aninhados que iteram sobre os índices das strings de entrada, enquanto que para o local há também dois loops aninhados, porém um para iterar sobre os índices da segunda string e outro para comparar as duas substrings. Para o exaustivo, pode-se afirmar que temos 3 loops aninhados, sendo entre eles um para cada indice de string e outro para comparação das subsmtrings geradas, ou seja, o algoritmo chega a ter complexidade O(n³).

Apesar de ter a mesma complexidade, o algoritmo local possui maior tempo de execução quando comparado ao Smith-Waterman. Isto pode ser explicado uma vez que o primeiro é executado repetidamente para aumentar as chances de encontrar um score melhor.

Alguns ajustes podem ser feitos no código exaustivo para que o tempo de execução seja menor. Um dos problemas é que para cada k (tamanho de substring) as substrings já analisadas são guardadas em um vetor. Isso para substrings pequenas são muito vantajosas, pois o as possíveis combinações da subsequência de DNA são poucas. Porém os scores altos são encontrados quando k possui valores mais altos e estes scores altos eliminam a necessidade de avaliar as combinações com k pequenos, ou seja, esta modificação que seria útil em um contexto onde k é pequeno, no final não é utilizada, mas as operações de guardar os elementos, fazer as comparações e limpar o vetor ainda são efetuados.

Com relação aos pontos de lentidão presente nos códigos, é possível afirmar que seriam os próprios loops são os agentes causadores. Isto será tratado a seguir com as técnicas de paralelização.

Paralelismo

Foi pedido que os alunos escolhessem entre os algoritmos anteriores para efetuar uma refatoração do código para que este fosse paralelizado. A paralelização seria feita em CPU utilizando a biblioteca OpenMP. Foi escolhido, então, o algoritmo de Busca Local esta tarefa, para que fosse possível explorar a geração de números aleatórios aprendidos em aula. Gerar números, de modo paralelo, aleatoriamente se torna um desafio em C++, pois há certos cuidados que devem ser tomados para que o programa não acabe gerando a mesma sequência de números ou que esta etapa não se torne um gargalo durante a execução.

Busca Local

A estratégia utilizada foi criar um gerador de números pseudo-aleatórios para cada thread:

std::vector<std::default_random_engine> engines;
std::random_device random;
std::default_random_engine generator1(random());
engines.push_back(generator1);
std::default_random_engine generator2(random());
engines.push_back(generator2);
std::default_random_engine generator3(random());
engines.push_back(generator3);
std::default_random_engine generator4(random());
engines.push_back(generator4);

Para gerar os números:

distribution(engines[round%4])

O ponto para efetuar o paralelismo foi no loop externo através de um "omp parallel for", onde calcula-se o tamanho das substrings. Além disso, foi criado também um ponto com "omp critical" para que não houvesse problemas de acesso ao tentar escrever na variável de score máximo.

Uma das iterações do codigo utilizava-se de funções para atividades especificas como gerar a substring, k ou p. Nesta iteração, foi observado que o tempo de execução era maior que o algoritmo sequencial. Quando os comandos executados pelas funções foram transferidos para dentro do código principal, houve uma melhoria substancial. A hipótese que tenta explicar o observado é que na chama de uma função, há um tempo adicional durante a execução devido a comandos extras relacionados a pilha de chamada, além de precisar alocar os argumentos nos endereços corretos. Nesta iteração, as funções criavam os geradores de números também, ou seja, para toda geração de índice, k, p ou substring uma nova variável de distribuição e um objeto gerador de números aleatórios era criada, aumentando consideravelmente o tempo de execução.

Abaixo estão alguns gráficos mostrando os resultados de tempo do algoritmo paralelo e outros fazendo a comparação com o sequencial. A quantidade de vezes que k era gerado foi o mesmo para ambos o sequencial e paralelo para que um comparativo fosse válido.

Foram feitos mais testes, porém estes possuem pares de DNA com tamanho muito maior (até 10000).

Como mostra os dois gráficos de comparativos de tempo, há uma melhoria de aproximadamente 40%. Uma suposição lógica, porém ingênua de se fazer, é assumir que o código seja 4 vezes mais rápido, pois está sendo executado em 4 threads ao invés de uma. Entretanto, deve-se considerar que além do comando "omp critical" transforma a área indicada em um algoritmo sequencial, tornando aquele pedaço em um gargalo durante a execução.

Busca Exaustiva

Para o paralelismo com GPU, foi pedido que o algoritmo exaustivo fosse paralelizado, porém, diferente do formato sequencial, a utilização da heurística de Smith-Waterman para o cálculo do score era obrigatoria. Isto acarretou em um comparativo inválido, mas que para este relatório foi desconsiderado. O algoritmo segue da mesma maneira do sequencial, porém o cálculo de score é feito com Smith-Waterman. Foi utilizado a biblioteca Thrust da NVIDIA. Vale ressaltar que a placa de vídeo utilizada foi dividida entre os vários alunos através de virtualização.

Houveram duas versões para a exaustiva paralelizada. A primeira consistia em paralelizar apenas o cálculo do score, enquanto a segunda tentava combinar a biblioteca OpenMP com a biblioteca Thrust. A segunda versão possuia o intuito de paralelizar o disparo de comandos para a gpu, através da utilização da OpenMP. Apesar de gerar resultados mais rápidos, foram encontradas inconsistências com relação ao score. Isso se deve ao fato de que o código criava apenas um vetor com o tamanho m. Apesar das execuções dos comandos Thrust serem sequenciais, ou seja, há apenas um comando sendo executado por vez, não havia garantia de que os comandos seriam enfileirados na sequência correta, podendo acarretar em influências do cálculo de um par de subsequência em outro. Por estes motivos, esta versão acabou por ser descartada para a análise em questão.

A outra versão que utiliza-se apenas da biblioteca Thrust tem como objetivo paralelizar o cálculo das linhas do Smith-Waterman. Isto é feito através da repetição de três etapas e uma outra etapa preparatoria:

  1. Calcular máximo entre diagonal e inserção.
  2. Calcular o máximo entre a etapa anterior e a deleção
  3. Obter valor máximo da linha para guardar o score.

Etapa preparatoria: Zerar os valores da linha inicial.

Foram criados 4 vetores na GPU: Duas para o score e duas para as strings. As strings de entrada foram transformadas em vetores de char e então enviadas à GPU. A etapa preparatória é feita sobre o vetor de score 1 antes do início de cálculo de um par de substrings.

O gráfico acima mostra uma inconsistência no tempo de acordo com o tamanho da entrada perto dos 60 caracteres. Isso pode ser explicado pelo Branch and Bound aplicado no algoritmo. Este fenômeno não está muito obvio no sequencial, pois este último foi testado mais vezes e portanto foi menos suscetível a possíveis paradas devido ao Branch and Bound. O algoritmo paralelo não foi testado mais vezes devido a demora do tempo de execução da versão atual.

Como é possível observar, algoritmo paralelo para tamanhos de string pequenos é muito inferior ao sequencial. Um fator que eleva o tempo de execução é a transferência dos dados para a GPU. O algoritmo em análise passou por uma iteração na qual vetores que seriam enviados a GPU eram criados toda vez em que o código criaria uma nova substring a ser analisada. Após executar alguns testes, este problema foi resolvido enviando os vetores no tamanho m no início do codigo e ao invés de enviar os vetores, alterar os índices que mostravam o início e fim da substring a ser analisada.

Apesar disto, ainda não é possível explicar o porquê do tempo de execução ser tão alto. Uma hipótese é que os comandos de criar tupla a partir dos vetores e criar iteradores seja muito demandante em questões de tempo de execução. Além disso, não se pode descartar a possibilidade do algoritmo ser mais rápido para entradas muito maiores. Para este relatório, testes maiores (na casa dos 10000 caracteres) estavam sendo executados, mas não foi possível terminá-los antes do prazo de entrega deste relatorio.